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INTRODUCTION 


The  program  NSUP  is  written  in  a  form  which  allows  data 
to  be  entered  interactively  on  a  DEC-10  computer  system  or  equiva¬ 
lent.  However  the  program  can  be  easily  adapted  to  batch  input 
and  to  other  computer  systems  by  altering  the  input  and  output 
statements.  Also,  the  velocity  history  is  defined  by  step  func¬ 
tions  in  each  of  the  six  degrees  of  freedom  to  aid  in  computing 
results  in  the  frequency  domain.  This  restriction  can  be  easily 
removed  and  arbitrary  velocity  histories  specified. 

Groups  of  Subroutines 

One  main  program  (MAIN)  and  19  subroutines  make  up  the 
program.  The  nineteen  subroutines  used  in  NSUP  can  be  arranged 
into  groups  according  to  their  functions  as  outlined  below. 

1.  Input  and  Initialization 

Much  of  the  input  and  initialization  is  done  in  the  main 
program  (MAIN).  In  addition  the  main  program  calls 
subroutines  EXP,  SINP,  MTN,  AFSTN,  and  AFSR  to  aid  in 
this  function. 

2.  Special  Function  Evaluation 

To  save  computer  time  the  exponential  and  trigometric 
functions  are  computed  from  prepared  tables.  Subroutines 
EXF  and  SINQ  perform  the  required  interpolation. 

3.  Computation  of  the  Source  Strength  Distribution  and  its 
Time  Derivative 

Subroutine  BLANC  computes  the  source  strength  and  its 
time  derivative  at  the  center  of  each  panel  from  the 
velocities  and  accelerations,  respectively,  at  the  panel 
centers  relative  to  the  local  ( free-surf ace  induced)  flow. 


4.  Free-Surface  Computations 

Subroutines  ACPTR,  PRFR,  and  CFSR  all  involve  free- 
surface  related  computations. 

5.  Body  Computations 

Subroutines  EBD,  POTB,  POTST,  SELF,  MAT JN ,  GE,  GO,  SOLID, 
and  PREP  all  aid  in  the  computing  body-related  matrices. 

The  main  program  and  each  of  the  nineteen  subroutines  are  described 
below. 


Program  Description 

MAIN:  Initially  this  program  checks  to  see  if  certain  large 

arrays  which  depend  only  the  body  geometry  have  been 
computed  and  stored  in  their  assigned  files  on  previous 
runs.  If  so  they  are  read  in  directly  and  do  not  have 
to  be  computed.  The  main  program  then  calls  subroutine 
EDB  which  reads  in  the  panel  description  of  the  hull  and 
computes  arrays  E  and  PX  if  they  are  not  available.  The 
matrix  E  gives  the  source  distribution  vector  (or  its 
time  derivative)  from  the  vector  specifying  the  normal  rela¬ 
tive  velocities  (or  accelerations)  at  the  panel  centers. 

Array  PX  gives  the  generalized  force  vector  induced  by 

3<J>BD 

a  source  distribution  vector  for  finite  speed  (the  pU—v. 

O  X 

term  in  the  pressure).  If  required  the  matrix  P,  giving 
the  force  vector  induced  by  the  time  derivatives  of  the 
source  strengths  of  the  panels,  is  computed  by  subroutine 
POTST.  Next  subroutine  MTN  is  called  to  specify  the 
body  velocities  for  impulsive  motion  in  each  of  the  six 
degrees  of  freedom.  Then  AFS  and  AFSR  subroutines  are 
called  to  define  and  initialize  the  f ree-surf ace.  After 
initializing  the  free-surface  the  main  program  calls  a 
pair  of  routines  to  compute  trigometric  and  exponential 
tables  (SINP  and  EXP)  for  later  use.  The  main  program 
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executes  the  time  stepping  loop  for  NTM  values  of  time. 
The  first  step  is  impulsive  and  covers  a  time  interval 
of  zero  duration  to  obtain  impulsive  pressures.  The 
forces  and  moments  computed  by  this  initial  time  step 
are  the  impulsive  values  resulting  from  starting  the 
body  with  finite  speed.  All  later  time  steps  are  of 
fixed  duration  (DT)  with  forces  and  moments  acting  on 
the  body  giving  a  time  history  over  the  duration  of  the 
calculation. 

The  first  subroutine  to  be  called  in  each  time  step  is 
ACPTR  which  computes  the  panel  pressures  and  normal  accelerations 
induced  by  the  free  surface.  Next  subroutine  ZBLANC  subtracts  the 
free-surface  component  of  normal  acceleration  as  computed  by 
ACPTR  from  the  accelerations  of  the  body  to  obtain  the  relative 
acceleration  at  each  panel  center.  It  then  applies  matrix  E  -to 
find  the  time  derivative  of  the  body  source  distribution. 

The  free-surface  induced  force  vector,  PF,  is  then 
computed  from  the  pressure  distribution  by  subroutine  PRFR  (called 
by  ACPTR).  Subroutine  POTB,  called  by  the  main  program,  computes 
the  force  vector  PB,  induced  by  the  time  rate  of  change  of  the 
panel  source  distribution  (through  matrix  P)  and  the  spatial 
derivative  of  the  body-induced  potential  in  the  x  direction  for 
finite  forward  speed  (through  matrix  PX).  The  two  force  vectors, 

PF  and  PB,  are  added  to  give  the  total  generalized  force  vector  PT 
for  each  of  the  six  degrees  of  freedom.  Finally  subroutine  CFSR 
is  used  to  advance  the  free  surface  by  a  single  time  step,  complet¬ 
ing  the  loop. 

Subroutine  EXP:  Sets  up  a  table  for  the  exponential  function. 

Subroutine  EXF:  Uses  the  table  to  compute  the  exponential  quickly. 

Subroutine  MIN:  Reads  in  impulsive  velocities  for  each  degree 
of  freedom,  time  step  size,  and  number  of  time  steps. 

Subroutine  SINP:  Sets  up  a  table  for  trigometric  functions. 

Subroutine  SINP:  Uses  the  table  to  compute  sine  and  cosine  func¬ 
tions  quickly. 
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Subroutine  ZBLACN:  Applies  the  specified  body  accelerations  for 

the  six  degrees  of  freedom  to  compute  the  resulting  nor¬ 
mal  accelerations  at  the  panel  centers.  The  free- 
surface  induced  normal  accelerations  ACNW(J)  are  sub¬ 
tracted  to  obtain  the  net  normal  acceleration  at  each 
panel  center,  ACN(J),  of  the  body  relative  to  the  fluid. 
These  accelerations  must  be  cancelled  by  the  time  deriv¬ 
ative  of  the  body  source  density  distribution,  ST(K). 

The  ACN(J)  vector  is  multiplied  by  the  E  matrix  to  get 
the  necessary  net  rate  of  change  of  the  panel  source 
densities,  vector  ST(K).  The  total  source  strength 
densities  are  accumulated  in  STOLD(K). 

Subroutine  ACPTR:  Computes  free-surface  induced  accelerations 
ACNff ( J )  and  pressures  PRFS(T)  at  panel  centers. 

Subroutine  PRFR:  Computes  generalized  force  vector  PF  for  six 

degrees  of  freedom  from  the  computed  free  surface  pres¬ 
sure  distribution  at  the  center  of  each  panel,  PRFS(T). 

Subroutine  AFSIN :  Reads  in  parameters  defining  the  free-surface 
representation,  then  calls  AFSR. 

Subroutine  AFSR:  Sets  up  and  initializes  the  free-surface 
representation. 

Subroutine  CFSR:  Advances  free-surface  by  one  time  increment. 

Moves  body  relative  to  the  free  surface.  Adds  the 
changes  in  free-surface  elevation  induced  by  the  body 
sources  acting  over  the  time  increment  to  the  free- 
surface  representation.  Second  order  effects  in  time 
are  included. 

Subroutine  EBD:  This  subroutine  reads  in  the  (x,  y,  z)  coordinates 
of  each  of  the  four  corner  points  into  a  set  of  arrays 
XPT(N),  YPT(N) ,  ZPT(N).  Panels  are  identified  by  a  set 
of  four  integers  giving  the  array  positions  of  the  four 
corner  points  of  each  panel.  Panel  areas,  normals  and 
center  point  coordinates  are  then  computed.  Finally, 


the  E  matrix  giving  the  source  time  derivative  distirbu- 
tion  for  a  set  of  prescribed  normal  accelerations  is 
computed.  The  inverse  of  E  is  computed  first  by  the 
subroutine  GE  which  gives  the  acceleration  induced  at 
any  panel  center  point,  J,  by  a  uniformly  distributed 
time  derivative  of  source  strength  density  of  unit  magni- 
tude  acting  over  any  surface  of  panel,  JL.  Subroutine 
MATIN  inverts  E  to  obtain  the  desired  form.  Simultane¬ 
ously,  the  matrix  PX  which  gives  the  x  component  of 
velocity  at  the  center  of  panel  J  induced  by  a  source 
strength  of  unit  magnitude  distributed  over  a  panel,  JL, 
is  computed. 

Subroutine  POTB:  This  subroutine  is  used  to  compute  the  general¬ 
ized  body-induced  force  vector,  PB,  for  all  six  degrees 

*  of  freedom  generated  by  a  known  source  strength  distribu¬ 
tion  and  its  time  derivative.  The  matrix  P(J,K)  is 
multiplied  by  a  vector,  ST(K),  representing  the  time 
derivative  of  source  densities  of  the  panels  to  obtain 

*  one  term  of  PB.  Similarly  the  term  proportional  to 

forward  speed  is  computed  from  the  matrix  PX  and  a 
vector  representing  the  accumulated  source  densities, 
STOLID(K) . 

» 

Subroutine  POTST:  Calculates  the  matrix  P(N,J)  giving  the  net 

force  or  moment  for  the  Jth  degree  of  freedom  induced  by 
a  unit  time  derivative  of  source  strength  over  panel  N. 
j  Fundamental  to  this  is  the  need  to  compute  the  potential 

integrated  over  each  panel  area  due  to  a  uniform  source 
density  over  every  other  panel.  For  panels  which  are 
far  apart  relative  to  their  dimensions  this  value  is, 

I  for  unit  source  density,  simply  proportional  to  the  pro¬ 

duct  of  their  areas  divided  by  the  distance  between 
centers.  The  method  used  here  is  to  divide  each  panel 
into  a  large  number  of  small  subpanels  and  then  calculate 
i  the  result  numerically,  adding  the  contributions  of 

each  subpanel  under  the  assumption  that  their  separa¬ 
tions  are  large  relative  to  their  dimensions. 
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Subroutine  SELF:  This  subroutine  is  called  by  POTST  to  compute 
diagonal  terms  in  the  P  matrix. 

Subroutine  MATJN:  The  matrix  inversion  routine  used  by  subroutine 
EBD  to  invert  matrix  E. 

Subroutine  GE:  A  subroutine  called  by  EBD  to  compute  the  elements 
of  matrix  E  prior  to  inversion.  It  computes  the  velocity 
(acceleration)  at  field  point  (XF,  YF,  ZF)  induced  by  a 
source  density  (time  rate  of  change  of  source  strength) 
of  value  unity  distributed  uniformly  over  panel  J. 

Subroutine  GO:  A  subroutine  called  by  GE. 

Subroutine  SOLID:  Also  called  by  subroutine  GE  to  compute  the  solid 
angle  of  a  panel  relative  to  the  field  point. 

Subroutine  PREP:  Prepares  all  panels  for  the  GE  subroutine.  It  is 
called  by  EDB  prior  to  using  GE. 

Sample  Input 

As  an  example  of  the  input  required  to  execute  NSUP 
interactively  from  a  remote  terminal,  a  sample  input  is  given 
below.  Asterisks  have  been  placed  at  the  beginning  of  all  lines 
output  by  the  program  to  distinguish  them  from  inputs  supplied  by 
the  user.  Comments  have  also  been  added.  They  are  distinguished 
by  placing  them  in  parentheses. 

♦TYPE  1  TO  READ  IN  PREVIOUSLY  CMPTD  BODY  ARRAYS 
1 

♦TYPE  1  TO  OUTPUT  FREE-SURFACE  ACCELERATIONS  AND  PRESSURES 
1 

(this  option  gives  a  printout  of  the  free-surface  induced  accel¬ 
erations  and  pressures  computed  at  the  panel  centers  for  every 
time  step) 

♦BODY  HAS  BEEN  DEFINED  WITH  60  PANELS 
(start  subroutine  MTN) 
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♦HOW  MANY  TIME  STEPS  AFTER  INITIAL  STEP 
150 

♦WHAT  IS  THE  TIME  STEP  INTERVAL  -  ZERO  FOR  FIRST  STEP 

0.10 

♦BODY  VELOCITIES 

♦VELOCITIES  ZERO  TIME  NEGATIVE  -  CONSTANT  TIME  POSITIVE 
(program  assumes  a  velocity  step  function  at  time  =  zero) 

♦WHAT  IS  SURGE  VELOCITY  FOR  POSITIVE  TIME 

0. 

♦WHAT  IS  SWAY  VELOCITY  FOR  POSITIVE  TIME 

0. 

♦WHAT  IS  HEAVE  VELOCITY  FOR  POSITIVE  TIME 

1.0 

♦WHAT  IS  ROLL  VELOCITY  FOR  POSITIVE  TIME 

0. 

♦WHAT  IS  PITCH  VELOCITY  FOR  POSITIVE  TIME 

0. 

♦WHAT  IS  YAW  VELOCITY  FOR  POSITIVE  TIME 

0. 

(end  subroutine  MTN  and  start  subroutine  AFSIN) 

♦INPUT  ACCELERATION  OF  GRAVITY  AND  FLUID  DENSITY 
1.0 
1.0 

♦WHAT  IS  FORWARD  SPEED? 

0.0 

♦INPUT  MAX  X  LENGTH  SCALE 

2.2 

♦INPUT  MIN  X  LENGTH  SCALE 

0.20 

♦INPUT  MAX  Y  SCALE 

2.2 

♦INPUT  MIN  Y  LENGTH  SCALE 

0.20 

♦INPUT  MAXIMUM  TIME  SCALE 

8.0 
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(The  time  scale  is  a  trade-off  between  accuracy  over  the  later 
time  steps  and  computational  time.  Generally  it  can  be  set  to 
a  value  somewhat  less  than  the  duration  of  the  simulation  with¬ 
out  introducing  significant  errors.) 

(END  OF  INPUT) 

Because  of  its  size  and  complexity  the  panels  defining  the  hull 
are  input  from  a  prepared  file  named  PANIN.DAT.  The  definitions 
of  the  variables  which  are  input  and  the  formats  by  which  they  are 
read  can  be  easily  seen  from  the  instructions  and  comments  at 
the  beginning  of  subroutine  EBD.  Integers  and  floating  point 
numbers  are  input  with  formats  415  and  3F10.0,  respectively.  The 
first  line  in  the  file  defines  NPT,  the  number  of  corner  points, 
and  NPAN,  the  number  of  panels: 

READ  (23,101)  NPT,  NPAN  . 

(The  file  PANIN.DAT  has  been  assigned  to  unit  23  previously.) 

Next  the  coordinates  of  the  corner  points  are  read  in  by: 

READ( 23,100) (XPT(N) , YPT(N) , ZPT(N) ,N=1 ,NPT)  . 

Finally  the  panels  are  defined  by  giving  the  four  corner  points 
(in  a  clockwise  sense  when  viewed  from  a  point  outside  of  the  body 
along  the  panel  normal).  These  four  points  are  identified  by 
integers  specifying  the  position  of  the  corner  point  in  the  corner 
point  array  previously  read  in: 

READ (23 , 101 ) ( KK(N , 1 ) , KK(N , 2 ) , KK(N , 3 ) , KK(N , 4 ) , N=1 , NPAN)  . 
Output 

The  program  starts  its  output  on  the  terminal  with  two 

columns  of  floating  point  numbers  giving  the  x  wave  numbers,  kxn> 

and  their  increments,  Akx  .  These  columns  are  repeated  for  the  y 

n 

wave  numbers,  kym,  and  their  increments,  Akym.  For  the  case 
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defined  by  the  sample  input,  the  x  and  y  wave  numbers  are  identical 
with  19  values  each  starting  at  0.01563  and  ending  at  5.33523. 

The  total  number  of  modes  is  19  x  19  =  361. 

During  each  time  step  information  is  typed  out.  As  an 
example  the  output  during  the  initial  time  step  of  the  computation 
resulting  from  the  sample  input  is  listed  below. 

IMPULSE  FORCE  AT  T=0+ 

(only  for  initial  step) 

TIME=0. 0000000 


FREE  SURFACE  INDUCED  FORCES - 


0.000000  0.000000 

0.000000  0.000000 

(initial  values) 


0.000000 

0.000000 


ACCELERATIONS 

0.00000000  0.00000000  0.00000000  0.00000000  0.0000000 0 

(repeated  over  12  lines  total) 


PRESSURES 

0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 

(Again  repeated  over  12  lines.  Initially,  the  free-surface 
disturbance  is  zero  and  does  not  induce  any  pressure  or 
accelerations  at  the  sixty  panel  centers) 


BODY  INDUCED  F0RCES= 


-0.000030 

0.006556 

-1.015025 

0.000039 

-0.000109 

-0.000124 

TOTAL  FORCE  COMPONENTS 

-0.000030 

0.006556 

-1.015025 

0.000039 

-0.000109 

-0.000124 

(END  OUTPUT  ON  TERMINAL) 


In  addition  to  the  output  on  the  terminal,  a  file  named  TOTAL.DAT 
is  created  containing  the  times  and  the  total  force  components 
computed  at  each  time  step.  Also  if  the  arrays  E,  P,  and  PX 
must  be  computed,  they  are  stored  on  permanent  files  for  future 
runs . 

Program  Listing 

A  listing  of  the  programs  is  given  over  the  next  few 

pages . 


C  MAIN  PROGRAM 


cohewi/evepaik  i20>  ,y?.vk  iho>  .zpaiu  ico) ,  arha<  20  ,stc  120 ) , 

.ACIK  ’-£0)  ,ACi IN!  120 >  .  AIK  120 ,3)  ,  E(  120)  120 , 0)  ,  ?RF S(  120)  .STOLE!  120) 

.  ,PM<  123,6) 

C0raiCIT/'FS/AI2(C0,33)  , SS(  30,20)  ,CC(  30,30)  , 

.  DICK  30 )  ,  DKY!  GO )  ,  AIC£(  CO )  ,  AI2Y!  30 ) 

coukon/bdo/npt!  iso)  ,vpt(  igq) ,zpt<  iso)  ,wrf!  iso  ,nrfr<  iso) 

.  ,ICa  130,4) 

COMMON/ A/NP AN,  IIPT,  CEE,  FJIO,  HEX,  NKY,  EVE,  DT,  TUI.  UFND 


DIMENSION  PF!6) ,P3(6) ,PT(6) 

COMPLEX  EYE 
EYE3 (0.0, 1.0) 

TYPE  4442 

4442  FORMAT!  ’  TYPE  1  TO  PE  AO  III  PREVIOUSLY  CMPTD  BODY  ARRAYS  ’  ) 

ACCEPT  400,  USED 
400  FORMAT!  13) 

C  COMPUTE  CODY  MATRICES  E  AND  ?X 
CALL  EBD(NSKB) 

TYPE  1111 

1111  FORMAT!  ’  TYPE  1  TO  OUTPT  PREE-STIRFCE  ACCELERAT 1 0113  .MID  PRESSURES’ ) 
ACCEPT  400 ,  tlTVP 
IF(  NSIC3.  NE.  1 )  CALL  POTST 
IF(  NSIC3.  NE.  1 )  GO  TO  7334 
C  READ  PRESSURE  MATRICES  INTO  CORE 

CP  Ell  (  UII!T=  1 ,  F  ILS3  *  PT3V’  ,  ACCESS3  ’  SEQIIT  ,  DEVICE3  ’  BSX:  *  ) 

OPErK  UWIT=2,  FILE3  ’PESV  ,  ACCESS3  ’SEQIN’  ,  DEVICE3  ’ DSIC:  ’ ) 

DO  544  K= 1 , 6 

READ!  1 )  (  PI  J ,  X)  ,  J=  1 ,  IIP  AN) 

544  READ!  2)  ( PX(  J,  IO  ,  J3  1 ,  II? AN) 

CLOSE!  UNIT3 1 ) 

CLOSE!  UIIIT=2) 

7334  NTM=0 

TYPE  16, RP AN 

16  FORMAT! ’  BODY  HAS  BEEN  DEFINED  WITH' , 13 ,* PANELS’ ) 

C  OPEN  FILE  TO  STORE  TOTAL  FORCE  COMPONENTS 

OPEN!  Utl IT=  19  ,  FILE3  ’  TOTAL ’  ,  ACCESS3  ’ SECOUT-  ,  DEVICE3  ’  EST:  ’ ) 

C  IMPULSIVE  BODY  VELOCITIES  (CONSTANT  AFTER  IMPULSIVE  FIRST  STEP) 

CALL  MTN(  NTH,  VMMT,  VYMT,  VZMT,  VRLMT,  VPIJT.  VY1AJT) 

C  INITIALIZE  FREE  SURFACE 

CALL  AF3IIK  T,  ECU,  SHE,  SGY,  SIIY) 

C  INITIALIZE  SPECIAL  FUNCTION  GENERATORS 
CALL  SINP 
CALL  EXP 

C  SET  TIME  TO  ZEP.O  AND  BEGIN  LOOP  IN  TIME! JTM) 

TIM=0. 00 
TYPE  2626 

2626  FORILAT! ’  IMPULSE  FORCE  AT  T30+’> 

DO  157  JTMTI3  1 ,  NTM 

JTM3 JTMM 

IF  (JTM.EQ.2)  TYPE  2627 

2627  FORMAT!’  F( TILDA)  IN  TIME  DOMAIN—  RESPONSE  TO  STEP  FUNCTION  ’) 
TYPE  17, TIM 

URITE! 19,2121)TIM 

17  FORMAT! //*  TIME3 ’ ,F13.6) 

C  COMPUTE  PANEL  PRESSURES  AND  NORMAL  ACCELERATIONS  FROM  7.3. 

CALL  ACPTR!  PF ,  NTVP) 

C  FIND  SOURCE  STRENGTHS  OF  PANELS 
IF!  JTM.NE. 1)  GO  TO  10 

C  FIRST  STEP  IS  IMPULSIVE  VELOCITIES  .ARE  FINAL  VALUES 
CALL  ZBLACII!  VXMT,  VYIIT,  VZMT,  VRLMT,  VPMT,  VYNMT,  JTM) 

C  INITIALIZE  SOURCE  STRENGTHS 


DO  4176  J= ! , NPAN 
4176  STCLDf J)=ST< J) 

GO  TO  20 

C  ZERO  ACCELERATIONS  AFTER  FIRST  IMPULSIVE  TIME  STEP 

C  FIRST  STEP  IS  IMPULSIVE-  VELOCITIES  REACH  FULL  VALUE  INSTANTLY 
10  CALL  ZBLACIKO.  ,0.  ,0.  ,0.  ,0.  ,0.  ,JTN) 

20  CONTINUE 

C  COMPUTE  FORCES  CONTRIBUTED  BY  BODY  SOURCE  DISTRIBUTION 
CALL  POTB(PB) 

C  ADD  TREE-SURFACE  AND  BODY  INDUCED  PRESSURES 
DO  76B  NGGN= 1 , 6 
PT(  IIGGN)  =PF(  NGGN)  +PB(  NGGN) 

768  CONTINUE 
TYPE  -*433 

4433  FORMAT! ’  TOTAL  FORCE  COMPONENTS’) 

TYPE  2121, PT(  1) , PTt  2) .PTI3) ,PT(4) ,PT(3) ,PT(6> 

WRITE!  19, 2 121) PTC  1)  ,?T(2)  ,?T(3)  ,PT(4)  ,PT(5)  ,PT(6) 

2121  FORMAT! 3I£,  3F 15 . 6) 

306  CONTINUE 

C  ADVANCE  FREE  SURFACE 

IF<  JTII.EQ.  1)  GO  TO  157 
CALL  CFSR(JTH) 

IF(  JTtl.GT.  1)  TIII=TIM+DT 
157  CONTINUE 

CLOSE!  UN I T= 19) 

STOP 

END 


SUBROUTINE  EXP 

C  PREPARE  TABLE  OF  EXPONENTIAL  FUNCTION 
C0M?I0N/EX/EX3T(  100) 

DEI=ENP(-0. 10) 

DET=EXP(-1.0) 

DA=  1 . 0 
J=0 

DO  974  M= l , 10 
DB=DA 
DA=DA*DET 
DO  974  N= 1. 10 
J=J+1 
D3S DBsDEI 
SN3TC J)=DD 
974  CONTI HUE 
RETURN 
END 


subroutine  rmu  ntm, vi ,  V2 .  V3 , V4 .  vs , V6) 

common/ A/ri?  an,  npt,  gee,  nno.  irxx.  nky.  eye,  dt.  tih.  ufwd 

COMPLEX  EYE 

C  HEADS  IN  CODY  VELOCITY  STEP  FNCTN  (RIGID  BODY  NOTION) 

C  STARTING  AT  T 1 131- ZERO 
TYPE  10 

ACCEPT  ICO. NTH 
NTH*  NTII+  1 
TYPE  000 
ACCEPT  200. DT 
TYPE  9 
TYPE  6 
TYPE  1 1 
ACCEPT  200. VI 
TYPE  12 
ACCEPT  200, V2 
TYPE  13 
ACCEPT  2C0.V0 
TYPE  14 
ACCEPT  200. V4 
TYPE  15 
ACCEPT  200. V5 
TYPE  16 
ACCEPT  200, V6 

10  FORMAT!  ’  BOW  ILAKY  TINE  STEPS  AFTER  INITIAL  STEP’) 

100  FORMAT!  15) 

200  FORMAT!  F 10. 0) 

300  FORMAT!  ’  WHAT  IS  TIME  STEP  INTERVAL  -ZERO  FOR  FIRST  STEP’) 

11  FORMAT!  ’  I 'EAT  IS  SURGE  VELOCITY  FOR  POSITIVE  THE’) 

12  FORMAT!  ’  VilAT  IS  SWAY  VELOCITY  FOR  POSITIVE  THE’) 

13  FORMAT!’  WHAT  IS  IIEAVE  VELOCITY  FOR  POSITIVE  TIMS’) 

14  FORMAT!  ’  K1IAT  IS  ROLL  VELOCITY  FOR  POSITIVE  TIMS’) 

15  FORILAT!  ’  Vfp.\T  IS  PITCH  VELOCITY  FOR  POSITIVE  TIMS’) 

16  FORMAT!’  VIE AT  IS  YAW  VELOCITY  FOR  POSITIVE  TIME’) 

9  FORMAT! ’  CCDY  VELOCITIES  ’) 

6  FORMAT!’  VELOCITIES  ZERO  TIME  NECATI VE-COtlST  TIME  POSITIVE’) 
RETURN 
EIID 


SUBROUTINE  SHIP 
COMMON/SNC/SRSTC  <02 )  ,  CST(  402) 
0P  =  <j .  233  1C3307/400 .  0 
F*0 . 00 

DO  13  N= 1.402 
SNST(R)=SIR(F> 

CST( N) »C0S(  F) 

Fa  F+DF 
15  CONTI HUE 
RETURN 
END 


J 


SCT3R0UT ITTE  SITHKSS.CS,  ARG) 

cofiKorf/SNc/srrsT(‘io2)  ,cstm©2) 

AE=ARG*0. 159 154943Q92 

WAE*AE 

AE=  AE-NAE 

I7< AE.LT.O.O) AS*  1 .CO+AE 

ND=AE*400.0+0.30 

NC=ND+1 

DAE*  ( AE-tlD*0.  0023) *3 .  141392634 

D.4F=DAII+DAI1 

SNl=SnST(NC) 

csi=cjT(r:c> 

511= SN  1+DA7*<  CS1-S31*DA3) 
CS3CS1-DA7:"(  SH1+CS  l-DAE) 

retoiui 

END 


subroutine  zblacni  acbx,  achy. albz. accrl. acbp, acbyw,  jttd 

C  CdPTTT.S  SOURCE  S’rR£NGTRS  SATISFYING  DODY  ROUND ARY  CFOTN 

i;oi'i:oN^Ei)/>a3/v;(  iao» .  v?an(  !20> , zpajk  ieo> ,  area!  ‘2o>  ,st(  120) , 

.  ACN!  120)  .ACNW!  120)  .AIK  120.C)  ,E(  120)  ,P<  120,S>  ,PIVFS(  120)  ,STOU)(  12«> 
.  .PX! 120,6) 


COMMON/ A/NP AM. nPT, GEE,  RHO.NXX.  NKY,  EYE,  DT.  TIN. UFVD 
COMPLEX  EYE 

OPEN!  UH!T=21 .FILE3  *E* .DEVICE* ’ DSK: * , ACCESS* ’SEQIN'  ) 

D  TYPE  6 

D  6  FORMAT!  *  START  ITERATION  FOR  TINE  DERIV  OF  3D?  SOIT.CE  STRENGTHS ’ ) 
DO  1COO  J  *  1 ,  NPAN 
X=XPAU( J) 

Y-  YPAM!  J) 

Z*ZPAN( J) 

AC  Xs  ACBX-f  AC3P"Z-  AC3YV* Y 
ACY* ACOY+  ACDY W«X- ACBPXtfZ 
ACZ=  AC3Z+ AC3RL* Y- AC3P*X 

C  DOT  PRODUCT  OF  ROY  ACEL+  3DY  NRML  FREE-SRFC  INDUCED  NRNL  ACEL 
ACM  J)  =  - ACNW!  .1)  +  ACX-AIl (  J ,  1)  +ACY*AN!  J ,  2)  ♦ACZ^A.I (  J ,  3) 

C  INTEGRATE  TIME  DERIVATIVE  FOR  TOTAL  SOURCE  STRENGTH 
C  AT  START  OF  TIME  STEP 

IF(  JTM.  GT.  2)  STOLD!  J)  =STOLD!  J)  +DTV-ST(  J ) 

ST( J) =0.00 
1800  CONTINUE 
D  TYPE  60 

D  60  FORMAT!  ’ NORMAL  ACCELERATIONS’) 

D  TYPE  15 . <  ACtH  N) , N* 1 , NPAN) 

D  13  FORMAT! 1X.5FX2.6) 

457  REWIND  21 

C  READ  IN  INVERSE  BODY  MATRIX 
DO  15C0  J*  1 ,  tlP  AM 
READ! 2 1 )  !  E!  J J ) .  JJ* 1 .NPAN) 

AC*  ACM  J) 

DO  14C0  XM.NPAN 
1400  ST!IO=STaO+E(IO=PAC 
15C0  CONTINUE 
RETURN 
END 


subroutine  acptri  pf  ,  ntpao 
c  computes  i.aye  ;ccelep.atio;:s  aiid  pressures  at 

C  PANEL  CENTERS 

CO.HION-'BD'XPAIN  !20)  ,  VP.uk  120)  ,  ZPAIK  120)  ,AREA<  !20)  ,3T(  120)  , 

.  ACIK  120)  .ACNVC  120)  ,  AIK  120,3)  ,  E(  120)  ,P(  120,0)  , ?UFS(  120) ,3T0LD<  120) 
. ,PX(  120,0) 

C0MK0tL'FS/AKZ<30,30> . SS( 30 , 30) . CC(  30,30) , 

.  DICK  30)  ,  DICY<  GO)  ,  AKXt  CO)  .  AXY<  30) 

C0I3I0N/rsi/A( 30, 30) ,8(30,30) ,AS< 30, GO) ,33(30, CO) 

COMMON/ A/NPAN .  HPT,  CEE, UNO,  NKX.  NXY,  EYE.  DT,  TIM,  U7WD 

COMPLEX  EYE .  SC ,  OX,  CX.  DY ,  CY,  3YC0H .  SCOII .  SX.  SKOH ,  3 1 ,  B2 ,  C 1 ,  C2 
COMPLEX  A. B, AS. OS 

DI  METIS  ION  DX(  30)  ,  CX(  30)  ,  PE(  6)  ,DX(30) 

DO  1300  J*  I ,  TIP  All 
AX*  AIK  J .  1 ) 

AY»AN(  J ,  2) 

AZ*AN(  J,3> 

X’XPAIK  J) 

Y*  YPAIH  J) 

Z* ZPAIK  J) 

ACT* 0.00 

PUT* O. 00 

DVDZ*0 . CO 

DO  16  II*  1.  NKX 

ccxx*  aek  rn  *x 

CALL  SIIIQ(S.C.CCXX) 

3X(  N)  *CI1PLX(  C.  S)  wDXX(  N) 

dik  ro  *  -  aick  id  *a:ck  X)  wax 
16  CX(  ri)  =  AKX(  to  v: AX- EYE 
'  DO  163  M*l,nKY 
CCYY*  AKY<  II)  wY 
C.VLL  S I  :iO<  S ,  C .  CCYY) 

EY*CIIPLX<  C.  S)  »DKY<  ID 
BYC0N=C0HJG(  BY) 

CY*  AKY(  II)  WAY- EYE 
DO  163  N=1,NXX 
ARGZ*  AKZ<  n .  ID  WZ 
CALL  EXF<  DEP , ARGZ) 

B I * DEPwflYwBJK  N) 

B2*  DEP* BYCON BX(  IO 

cz*akz(H.M)waz 

SC*BlwA(N.It> 

SC0N*B2*AS(  N,  ID 
SX*31wB<  II, II) 

S:C0N*B2wB3(  II.  M) 
c  i  *cx<  :o  +cY-cz 

C2*CX(rt)-CY-CZ 
ACT*  ACT-C 1 "SC-C2WSC0N 

I; ;  X*  D VZ+EV.D ;AICK  10  *(  ClwSX<-C2*SK0N)  /SQRT(  .VKZ(  N ,  M)  ) 

PRT*  PRT+SC+SCON 
163  CONTINUE 

C  NOIUIVL  ACCELERATION  INDUCED  BY  FREE  SURFACE 
OVDZ*  DVDZw'JFI.T)w32RT(  CEE)  wDVDZ 
\CNV(  j) *  ACTwCEE+DVBZ 
C  PRESSURE  INDUCED  NY  FEZ Z  SURFACE 
PES'  .J)  sPRTwGEEwFZIO 
1300  CONTINUE 

C.VLL  PRFR(PF) 

CO  rONI'ATI  •  ACCELERATIONS’  ) 

60  FCRMATf  *  PRESSURES ’ ) 

1070  FORMAT' 1X.3F16.Q) 

i "(  NTPAC . HE . 1 )  RETURN 


TYPE  B* 

tyre  ’  070 , (  \CNV(  J)  .  J*  1 .  NP  VN) 

TV  ^  ^0 

TYPE  '070,  t  PI1FS(  J) ,  J*  1  ,NP.\N) 

RETURN 

END 


r 


SUBROUTINE  PRFIUPF) 

C  COMPUTE  FORCES  Ai!D  nOiEHTS  DUE  TO  FREE  SURFACE  INDUCED  PRESSURES 
COrUION/BD/XPAIK  120)  ,Y?AH<  120)  ,2?Ali(  120)  , AREA!  :20)  ,ST(  120)  , 

.Acrr<  120)  ,Acrm<  120)  ,aik  120,3)  ,e<  120)  ,p<  120,0  .prfsc  120)  .stoldc  120) 
.  ,PXC 120,6) 

COMMON/ A/NP AN ,  HPT,  GEE ,  IUIO ,  NICF.  NKY ,  EYE ,  DT ,  T I M.  UFND 

DIMENSION  PF( 6) 

COMPLEX  EYE 

Xl=0.0 

X2=0.GG 

X3=0.0 

X4=  0 . 00 

X5=0.CO 

X6=Q . 00 

DO  723  J= 1 . HP AH 
?RS*PRFS( J) wAREAC J) 

FPJ{=-AN(  J,  l )  *PE-3 
FRY=-AII(  J,2)*PRS 
FRZ*-.'\IH  J,3)*PIU3 
XF*XPAH( J) 

YF=  YP  AIK  J) 

ZF=ZPAN( J) 

Xl=Xl+FRX 

X2=X2+FRY 

X3=X3+FRZ 

N4=  X4+YF*FRZ-ZF*FTVY 
N5=  X3+ZF«FRX-XF-::FP,Z 
X6=X6+XF*FRY-YF«FRX 
723  CONTINUE 
TYPE  80 

80  FORMAT! ’  FREE  SURFACE  INDUCED  FORCES - ’ ) 

TYPE  <  0 ,  XI ,  X2 ,  X3 ,  X4 ,  IS .  X6 

PF( 1)=X1 

PF(2)=X2 

PFC  3) =X3 

PF( 4) =  N4 

PFI  5)  =XS 

PF<  6) =X6 

40  FORMAT!  5X,  3F 15 . 7) 

RETURN 

END 


SUBROUTINE  AFS IN(  T.30X, SXX,  3GY,  SKY) 
c  input  free  surface  parameters 

COKTCOIL'FS.'AKZI 30 . 20)  ,SS(30.30)  , CCC30.C0)  , 

.  DICK  30)  ,  DKY<  30)  ,  AICK  GO)  ,  AICY(  CO) 

comkon/a/npan ,  ijpt ,  gee ,  rho ,  nice  nicy,  eye ,  dt,  tim,  tjfwd 

COMPLEX  EYE 
TYPE  171 

ACCEPT  1 00 .  GEE ,  PJJO 
TYPE  200 
ACCEPT  100,  UFVD 
TYPE  180 
ACCEPT  100, E GX 
TYPE  181 
ACCEPT  100, S MX 
TYPE  132 
ACCEPT  100. EGY 
TYPE  ISO 
ACCEPT  100, SKY 
TYPE  20 
ACCEPT  100, T 

CALL  AFSR( T, BGII, 3KX, BGY,  SKY) 
type  20 1 ,  ruci.  NICY 
201  FORKAT(  211,  215) 

TYPE  21 

BO  132  N=  1 .  RICH 

152  TYPE  102 ,  AICK ID  ,  DICK  N) 

TYPE  22 

DO  153  HM.rilCY 

153  TYPE  102 ,  AICYC  ID  ,  DXY(  fl) 

171  FORMAT!  '  INPUT  ACCELERATION  OF  GRAVITY  AND  FLUID  DENSITY’) 
100  FORMAT!  F 10.0) 

180  FORMAT!’  INPUT  MAX  X  LENGTH  SCALE’) 

131  FCPulAT!  ’  INPUT  KIN  X  LENGTH  SCALE’) 

132  FORILAT!  ’  INPUT  MAX  Y  SCALE’) 

-  ISO  FORMAT!  ’  INPUT  Kill  Y  LENGTH  SCALE’) 

20  FORMAT!’  INPUT  MAXIMUM  TIKE  SCALE’) 

21  FORMAT!’  X  NAVE  NUM3EP.S  ’ ) 

22  FORMAT!’  Y  NAVE  IIUII3EIIS  ’ ) 

102  FORMAT!  1X.5F12.5) 

200  FORMAT!  ’  WHAT  IS  FORWARD  SPEED?’) 

RETURN 

END 


~-:t.outi:;s  t,  sgx,  s:n,  sgy,  stiy) 

C  IIITI IALI2E3  FREE  SURFACE 

Cf.ilKCtl/BD/lPAIiC  120)  ,  7? AIK  120)  ,  7-?AH(  120)  ,  AREA!  120)  ,3T(  120)  , 

. i\cnc  iso?  ,Acnx  :so)  ,aiu  120,3)  ,21 120) ,?< 120,0)  ,?n73<  1201 .stcldi  120) 

. ,?XC 120,6) 

C0MKOn/FS/'AKZ(3Q,30) ,SS(30,30> ,CC(3Q,GQ)  , 

.  DICK  30)  ,  DICY(  GO)  ,  Area  GO)  ,  AICY (  GO) 
cern:o;jxFsi/A(CO,co)  .bigo.go)  ,  as< 30, go)  ,3S(30,oo) 

CQIRIOW/'A/'HPAH ,  HPT ,  GEE ,  RHO ,  KXX,  RKY,  EYE,  DT,  Till,  UF’VD 

COMPLEX  A,  B. AS, ES, EYE 
DKTT*  1 . 0/-(  CEESTST) 

ITBX*  0 . 5/  (  BGXsDKTT) 

RBY*  O.S/(  BGY-DXTT) 

AIK*  1 . 0/SMX 

amy=  1 .  o/siiy 

rr=o 

11=9 

ACLD*  0 . 00 

1  IF(R.GE.RBX)  GO  TO  2 

rr=  rr+ 1 

AKX(IO=N*H*DICrT 

A0LD=AK3i(N) 

IF(  AOLD.GE.AIK)  RBX*R 
GO  TO  1 

2  r=r+i 

AKX<  N)  = AOLD+ 1 . 0/3GX 
AOLD=AKX(N) 

IF( AOLD.LE. AIK)  GO  TO  2 

nicx=r 

Aina  rrxx+ 1 )  =  aold+  i  .  o/bgx 

AOLD=  0 . 00 

11  IF(H.GE.IIBY)  GO  TO  22 
II=M+1 

AKY<  ID  =  M*It*DXTT 
AOLD=AKY(ID 
!F(AOLD.GE.AIIY)  RBY=M 
GO  TO  11 
22  II*  11+ 1 

A1CY(  ID  =  AOLD+ 1 . 0/EGY 
AOLD* AKY<  ID 

IF < AOLD.LE. AHY)  GO  TO  22 
mcY=n 

AXY<  M+  1 )  =  AICY<  M)  + 1 . Q/EGY 
AFF=  0.000 

no  i7  r=i,rkx 

AF=  t  AKX(  N+  1 )  +AKXC  H)  )  *0 . 50 
DKX<  tl)  =AF-AFF 

17  AFF* AF 
BFF=0.OO 

BO  18  Rsi.RKY 

BE*  (  AKY<  rt+  1 )  + AKY(  ro  )  *0 . 50 

B!CY(R)*EF-3FF 

18  OFF*  BF 

DO  ICO  R*  1 ,  rkx 
do  ioo  n*i,rucY 
Adt.IDMO.  0.0.0) 

3c:r.rD*<o. o.o.o) 

A3<R.ID*<0.0.0.0) 

cs<;<,ii)*(o. o.o.o) 

AXZ<  H ,  TD  *  SORT!  AICX1  IT)  s*2+ AICYC  M)  **2) 

S I G*  SCim  GEE* AKZ (  R ,  ID  ) 

SS<  R,  ID  *S IR( S IG*DT) 


CC( R,  ID  »COS( S IG*DT) 
100  CORTIROE 
R2TURR 
ERD 


SUBROUTINE  CFSRC  JTM) 

C  ADVANCES  FREE  SURFACE  WAVE  SPECTRA  IN  TIKE 

COMMON/ BD/XPANC  120)  ,  VP  AIK  120)  ,  ZPANC  :20)  ,AR£AC  120)  ,ST(  120)  , 

•  ABN (  120)  ,  ACINIC  120)  ,  All C  120,0)  ,  E(  120)  ,  ?(  120,6)  ,?RFS<  120)  ,  STOLDC  120) 
PIN  120,6) 

C0f!K0rr/FS/AKZ(30,S0)  ,  SS(  30,  CO)  ,  CC(  30 , 20)  , 

.  TICK  30)  .  niC/(  30)  ,  Aide  30)  ,  AIC/C  30) 

COKCION/FS  1/AC 30,30)  ,B(30,30)  , ASC 30, 30)  ,BS(00,30) 

COKIIOII/ A/NP AH ,  NPT,  GEE,  R1I0,  NKX,  NKY,  EYE,  DT,  Till,  UFTvD 

01  KENS  ION  CXC  30) 

COUP LEX  EYE , AT , ATS , CX, CY, CXY1 , CXYS ,  A.  B ,  AS ,  BS , DFUD 
C  FIRST  TIME  STEP  IS  PURS  IMPULSE-  TIHE=0.0  AFTER  FIRST  TIME  STEP 
IFC JTH.EO. 1)  CO  TO  6 
DO  100  N=1,NKX 
CXX=  AKXC  N)  -•kUFUD-'.’.-DT 
CALL  S I RO(  S , C , CXX) 

DFVD=CHPLXCC,-S) 

DO  100  11=1,  NIC/ 

CT=CC<  IT,  II) 

STT=SS( N, K) 

AT®  A(  N ,  ID  ::=CT+B(  II ,  M>  *STT 
ATS®  ASC  ri ,  II)  :::CT+  ESC  If,  II)  i':STT 
PCN,M)=BCN,M)  =::CT- AC  N ,  II)  *STT 
BSC  N ,  ID  =3S(  N .  ID  ::-CT-AS(  N ,  M)  =:<STT 
A<  N ,  K)  =  AT 
AS  ( N ,  ID  =  ATS 

C  HOI'S  FREE  SURFACE  RELATIVE  TO  BODY  WITH  FIN)  SPEED 
AC  II,  II)  -  AC  N,  II)  ADFND 
BC  N ,  ID  =  DC  IT ,  M)  <;DFUD 
AS(  N,  II)  =  ASC  N,  ID  sBFT/D 
BSC  N,  II)  =  BS  C  N ,  II)  :;:DFTID 
100  CONTINUE 
6  CONTINUE 

C  NOW  ADD  EFFECTS  OF  SOURCE  PANELS  ACTING  OVER  ONE  TIME  STEP 
CO  1500  J®  1 ,  RPAII 

STAR®  ( 3TC  J)  "0 . 5 0;." 0T+  STOLD (  J) )  *DT 
C  ST  IS  TIME  RATE  OF  CRINGE  OF  SOURCE  STRENGTH 
C  STOLD  IS  SOURCE  STRENGTH  AT  START  OF  TIME  STEP 
C  STAR  IS  AVERAGE  VALUE  OF  STRENGTH  OVER  TIME  STEP 
STAR2® STOLDC J) »DT 
IFC  JTII.EQ.  1)  STAR.=  ST(  J) 

IFC  JTII.EQ.  1 )  STAR2=0 . 00 
STAR=  STARSt AREAC  J)  SO . 6366  197724 
STAR2=STAR2*AREAC  J)  S© . 3 133 1 
X=XPAHCJ) 

IFC  JTH. GT. 1 )  X=X-UFVD*DT* . 50 
Y=YPANC J) 

Z= ZPANC  J) 

DO  93  N= 1 , NKX 
CXX®  AKXC  N) «X 
CALL  SINQCS,C,CXX) 

93  CXC  N) =CMPLX(  C,  -S) 

DO  94  11=1,  NICY 
CYY=AICY(  IDa-Y 
CALL  S INCH S, C, CYY) 

C Y®  CMPLXC  C , -S ) 

CO  94  N® 1 . NKX 
ARGZ®  AKZC  N.IDSZ 
CALL  EXFCDEP, AROZ) 

CXYI  ®  CXC  N)  s:CY*DS? 

CXYS® CXC  ID  sDEPsCQNJGC  CY) 

AC  N ,  ID  =  AC  II ,  ID  +STAJV:CXY1 
ASC  N ,  II)  =  ASC  N ,  ID  +STAR£CXYS 

C  BC N -  II)  AND  ESC H,  ID  INCREMENTED  NEGLECTINC  CHANGES  IN  SOURCE  STRENGTH 


C  OVER  TIKE  INTERVAL  APE  SECOND  ORDER  IN  DT 

B(  rr,  id  =  bc  ri ,  in  -starsscxyi  *ssc  n.  id 

BSC  N ,  ID  =  33t  N,  ID  -STARHsCXYSsSSC  N ,  H> 

94  CONTINUE 
1500  CONTINUE 
RETURN 
END 


SUBROUTINE  EBDC  IISIG3) 

C  IHTITIALIZE  PANELS  AND  COMPUTE  BODY  MATRIX 
C 

CONMON/DD/XPAITC  120)  ,  YPAIK  120)  ,ZPAN<  120)  ,  AREA(  120)  ,ST(  120)  , 

.  ACIK  120)  .ACIR/C  120)  ,  AIK  120,0)  ,  EC  120)  ,  PC  120,6)  ,  PRFSC  120)  ,STOLD<  120) 
.  ,?I1C  120,6) 


C0MM0N/BD2/XPTC  130)  ,  YPTC  130)  ,  ZPTC  150)  ,1/RFC  150)  ,  WBFRC  130) 

.  ,ICCC  150,4) 

COMMON/ A/NP AN ,  I1PT,  GEE,  RIO,  NEC.  NKY,  EYE,  DT,  TIM,  UFI/D 
DIMENSION  EP< 120) ,E??( 120) 

COMPLEX  EYE 

C  READ  IN  BODY  PANEL  PARAMETERS 

OPENC UNIT=20 , F  ILS= ’ PAN ! N’ ,  DEVICE*  ’  BSICs 1 , ACCESS3 ’ SERIN’ ) 

1C  I  FORHATC  4 15) 

ICO  FORMAT!  3F 1 0 . 0 ) 

C  NUMBER  OF  POINTS  AND  PANELS 
P-EADC  23.101)  HPT,  IIP  AN 
C  COORDINATES  OF  POINTS 

READC  23 .  100)  C  XPTC  ID  .  YPTC  N>  ,  ZPTC  H)  ,  N=  1 ,  N?T) 

C  DEFINE  CORNER  POINTS  OF  EACH  PANEL 

READC23,  JOIMXXOJ.  1),:CCCII,2)  ,  XICC  N,  3)  ,  XHC  N,  4)  ,11*1,  If?  AN) 

CLOSEC  UNIT3 20) 

C  COMPUTE  PANEL  AREAS 
DO  150  J= 1 , NP AN 
Kl  =  ia«  J,  1) 

IC2=KIC<  J.2) 

ICS  =  XICC  J,0> 

K4=KKC J.4) 

IFCK4.EQ.0)  GO  TO  3 

XPAIIC  J)  =  i  XPTC  XI )  +XPTC  !C2>+:CPT(  ICS)  +XPTI  X4>  )  «0 . 25 
'■PANC  J )  -  (  YPTC  XI )  +  YPTC  X2 )  +  YPTC  K3 )  +  YPTC  X4  >  )  WO .  25 
2PANC  J)  *  C  ZPTC  KI )  +ZPTC IC2)  +ZPTC  ICS)  +ZPTC IC4)  )  WO .  25 
GO  TO  9 

C  TRIANGULAR  PANELS 

8  I IP  AN C  ,!)=<.  XPTC IC1 )  +ICPTC IC2)  +  XPTC IC3)  1/3.00 
YPAHC  J)  =  C  YPTC  XI )  +YPTC IC2)  +YPTC IC3) )  /3 . 00 
ZPANC  J)  =  C  ZPTC  Kl )  +ZPTC  X2)  +ZPTC  X3) )  /3 . 00 

-  j’/g 

9  XA=  XPTC  K3)- XPTC  Xl) 

X3=  XPTC  IC4)  -XPTC  IC2) 

YA=  YPTCK3 )  -  YPTC  Id ) 

Y3=  YPTC  IC4 )  -  YPTC  IC2 ) 

ZA=ZPT(  ICO)  -ZPTC  XI ) 

Z3=  ZPTC  K4 )  -ZPTC  :C2) 

C  COMPUTE  PANEL  AREAS 
AZ=  XAW  YE- YAWX3 
AX- YAWZ3-ZAWYB 
A7=ZA=::X3-XAWZB 
ARE* SC1RTC  AX*AX+AY*AY+AZ*AZ) 

AREAC  J)  =AKEWO. 50 
.MIC  J,  1)=-AX/ARE 
ARC  J,  2)  =-AY/ARE 
AIIC  J,3)=-AZ/ARE 
150  CONTINUE 

806  FORILATC  ’  J '  ,  9X,  ’NX’  ,9X,  ’NY’  ,9X,  ’NZ’  ,9X,  ’X?’  ,9X,  ’  YP’  ,9X,  ’ZP"  , 

.  9X, ’AREA’ ) 

807  FORILATC  15 , 7F 1 1 . 4) 

I F  C  NSICB .  Ell.  1 )  RETURN 

C  RETURN  IF  E  AND  PX  ARRAYS  RIVE  BEEN  COMPUTED  IN  EARLIER  RUNS  T/IT3 
C  SAME  BODY 


OPEN!  UNIT=2 1 ,  FILE*  '  2’  , ACCESS* * S2G0UT'  ,  DEVICE*  ’  DS!C:  ’> 

<:?sm  ut: ir=2<- ,  fils*  *  ?:sv  ,  access*  •  segout’ ,  device*  •  qs:o  • ) 
jo  uoa  j=i.npan 
jjj*j 

CALL  PREPCJJJ) 

ST( J) =0 . 00 
DO  1308  K= 1 , 6 
PX(  J,  10=0.00 
1308  coirrirruE 

DO  308  J=  1 ,  NPAH 

JJJ=J 

AX*  AN! J,  1) 

AY*  AN!  J.2) 

AZ=AN! J,3) 

XF=XPAITl  J) 

VF=YPMI(  J) 

ZF=ZPAN! J) 

DO  157  L=l.riPT 

WRF!  L)  =SCR7(  (  KPT!  L)  -XF) **2+!  YPTC  L)  -Y7)  **2+(  Z?7(  L)  -Z7)  **2) 
157  KRFIU I.)  =SGET!  (  MPT!  L)  -XF) *« 2+<  YPT!  L)  -  VF)  **2+<  Z?T<  L)  +ZF>**2> 
DO  300  JL=  1 ,  IIP  AN 
JLJ= JL 

CALL  CE(  XT' ,  YF ,  ZF ,  JLJ ,  VX,  VY.  VZ,  VXR,  VYR,  VZR,  JJ J) 

VX=VX+VXR 
VY=  VY+VYR 
VZ* VZ+VZR 

C  COMPUTE  NORMAL  VELOCITY  AT  PANEL  J  DUE  TO  PANEL  JL 
E(  JL)  =  AX*VX+AY*VY+AZ*VZ 
C  INCREMENTT  PX  MATRIX 

Fill  =  - AREA (  J)  tfVXXAH!  J ,  1 ) 

F.T2»-AHEA(  J)  *VXSAN<  J,2) 

7  33  =  -  AREA'.  J )  *  VX*  AN  (  J .  3 ) 

PX! JL, 1 )  =  PX!  JL , 1)+FR1 
PX( JL, 2) =  PX!  JL, 2) +  732 
PX(  JL,  3)  =PX!  JL, 3) +FR3 
PX( JL, 4) *PM( JL, 4) +YF*FR3-ZF«FR2 
?"( JL, 3) =  PX! JL, 5) +ZF*FRl-XF»FR3 
?X(  JL , 6 ) =  PX(  JL . 6 ) +XF«FB2- YF*FR1 
309  CONTINUE 

T.TllTEt  2 1 )  ( E(  JL)  ,  JL=  1 ,  HP  AH) 

308  CONTIIiUE 

DO  2424  K=  1 , 6 

2424  WRITE!  24)  (  PX(  JL,  K)  ,  JL=  1 ,  NPAH) 

CLOSE! UN IT=21) 

CLOSE! UN IT=24) 

C  INVERT  E  MATRIX 

CALL  MATIN!  NP AN) 

RETURN 

END 


SUBROUTINE  POTB(PP) 

C  FIND  FORCES  AND  MOMENTS  INDUCED  3Y  THE  RATS  OF  CHANGE  OF  SRCE  STRNGTH 
C  IN  SPACE  FIXED  COORDINATES 

CCMHON/BD/XPAN!  120)  ,  VP  AIT C  120)  ,ZPAN(  120)  ,  AREA!  120)  ,ST(  120)  , 

. ACIH  120) , ACIW! 120) , AIK  120,3) ,E(  120) ,P( 120,6) , PRFS(  220) ,STOLD< 120) 

. ,PXC 120,6) 


COMMON/ A/ NP AN ,  NPT,  GEE,  RHO, NXX,  NKY.  EYE ,  DT,  TIM,  UFVD 

COMPLEX  EYE 
DIMENSION  PP( 6) 

PP( 1)=0.00 

PP! 2) =0 . 00 

PP( 3) =0. 0 

FP( 4) =0.00 

PP( 3) =0.0 

PP( 6) =0. 0 

DO  1500  J= 1 , NPAN 

C  ST( J)  IS  TIME  RATE  OF  CHANGE  IN  HULL-FIXED  SYSTEM  OF  SOURCE  STRENGTH 
C  OF  PANEL  J 

C  STAY  IS  AVERAGE  SOURCE  STRENGTH  OVER  TIME  STEP  AT  CENTER  OF  PANEL  J 
3TAV=  STOLD ( J )  +0 .  SsDT-sST!  J) 

C  FIRST  TIME  STEP  IS  PURE  IMPULSE 
IF! JTM.EG. 1)  STAV=  0 . 0 
C  TIME  DERIVATIVE  IN  SPACE  FIXED  SYSTEM 
DO  1200  K= 1 , 6 

PP(  X) =PP(  IO  +(  ST!  J) SP!  J ,  IO  -STAV*UFWD*PX<  J , K) ) *REO 
1200  CONTINUE 
1300  CONTINUE 
TYPE  1870 

1870  FORMAT! ’  BODY  INDUCE  FORCES  =  * ) 

TYPE  2020,  PP( 1) ,PP<2) ,PP(3) ,PP(4) ,PP(5) ,PP(6) 

2020  FORMAT! 3X, 3F 13 . 6) 

RETUI1N 

END 


SUBROUTINE  POTST 


i:0'~:o':/BD/:’?A;r(  ico  ,y?an<  120  ,zpan<  ico)  ,ap.sa<  !20>  ,-3T<  120) , 

. iiciK  120)  ,Ac:n;<  120) ,  aik  ico,3)  ,e<  120)  ,?<;  120,0)  ,?r.?3<  :20>  ,stold<  120) 
.  , PX( 120,6) 


coki:d:t/bd2/xpt<  iso)  ,ypt<  130)  ,zptc  iso  ,wrf(  150  .urfiu  iso) 

.  ,ICK(  130,4) 

COMMON/ A/NP AN ,  NPT.GEE ,  RHO ,  USX,  NKY,  EYE ,  DT,  Till,  UP’.® 
COMPLEX  A, B. EYE 

DIMENSION  XPSL( 3, 4)  ,XPSLR(3,4)  .  PBB(  120,  120) 
C0MK0N/PTST/ARE4C  200 , 4)  ,X4(200,4)  ,Y4(200,4)  ,Z<-(200,4) 

. , SEL( 200,4) 

BO  1500  J=  1 ,  nPAII 
ARE4( J,4)=-1.Q 
JT=4 

IF(ICK( J.4) .EQ.O)  JT-3 
CO  1500  JJ= 1 , JT 
J2=  1 

IF(JJ.LT.JT)  J2= JJ+ 1 
ICF=KK(  J,JJ) 

KG=KX(  J,  J2) 

X4(  J ,  JJ)  =  (  XPT(  Iff)  +XPTC  ICG)  + XPAIK  J) ) /3. 0 
Y4(  J ,  JJ)  =  (  YPT(  Iff)  +YPT(  ICG)  +YPAN<  J) )/3.Q 
Z4<  J , JJ)  =  <  ZPT<  Iff)  +ZPTC  ICG)  +ZPAN(  J) )/3.0 
AF=XPT<  Iff)  -XPAIK  J) 

BF=YPT(  Iff) -YPAIK J) 

CF^ZPTC  Iff)  -ZPAIK  J) 

AG=  XPT(  ICG)  -XPAIK  J) 

BG=  YPT(  KG)  -YPAIK  J) 

CG=  ZPT(  ICG)  -ZPAIK  J) 

CALL  SELF<  AF , BF , CF , AG, EG, CG, FEE) 

SEL( J, JJ)=FEE 
CR=.4F«BG-BF*AC 
AR=BF*CG-CF«8G 
EE=CF«AG-AF*CG 

ARS4(  J ,  JJ)  =0.3«SQim  AR«AR+BR*BR+CR*CR) 

1500  CONTINUE 

DO  127  HJ=  1 ,  KPAIJ 
DO  1277  MJ=1,NPAH 
1277  PB3(NJ,MJ)=0.00 
P(NJ, 1 ) =0.00 
P(NJ,2)=0.00 
P( NJ,3) =0.00 
P(  NJ, 4) =0.00 
P(  NJ, 5) =Q . GO 
P( NJ, 6) =0 . 00 
DO  123  NIC=1,4 
ARN=  ARE4(  NJ ,  NIO 
IF<  ARIT.  LT.  0 . 0)  GO  TO  128 
P1=0.00 
P2=0. CO 
P3=0. 00 
P4=0. 00 
P5=0.00 
P6=0.C0 
X-X4(  NJ ,  NIO 
Y=Y4(  iu.ino 
Z=Z4(  NJ ,  NIO 
DO  133  MJ=l,nPAIT 
BO  130  rEC=l,4 


YF=Y<-<MJ,MK> 

ZF3  Z4<  MJ ,  NX) 

ARM3  ARE  4  (  ru ,  fCO 

IF< ARM. LT. O.OO)  GO  TO  1G8 

IF(KJ.HE.MJ)  GO  TO  140 

IFIMK.NE.HIO  GO  TO  140 

FRA=SEL<  MJ ,  IEO  /AMI 

GO  TO  1330 

140  RA=SQRT(  ( X->IF)  SS2+(  Y-YF)  **2+( Z-ZF)  **2) 

FRA3  ARM/ RA 

1330  R3=SQRT<  <M-MF)s»2+(Y-YF)**2+(Z+ZF)**2) 

FRA3  FRA- ARM/RB 
FRX3-AH(MJ,  l)»FRA 
FRY3- AIT  (  MJ ,  2)  :.'FRA 
Fnz=-i\nau,3)aFRA 
Pl  =  Pl+FPJi 
P2=P2+FRY 
P3=P3+FRZ 

P4=  P4+YF*FRZ-ZF*FRY 
PS= P5+ZF*FRX-KF*FRZ: 

P6=P6+iaF*FRY-YF*FRX 
PBB(  NJ ,  MJ)  3  PBB(  HJ ,  MJ)  +FRA*ARH 
133  CONTINUE 

P(  ru,  i)=p<ru,  d+pi*arn 

PINJ.  2)3P(tU,2)-t-P2*ARH 
P(  RJ ,  3)  =P<  ru  ,  G)  +P3:.“-ARN 
P<  NJ ,  4)  3P(  ITJ,  4)  +P4SARN 
P<  riJ ,  3 )  3  P<  IIJ ,  3 )  +P5*ARN 
P(  riJ ,  6)  =  P(  ru ,  6 )  +P6SARH 
128  CONTINUE 
127  CONTI NUE 

OPENf  UNIT3 1 , FILE3 ’ PTSV’ , ACCESS3 ’ SEGOUT’ , DEVICE3  * DSX: • ) 
DO  534  KG= 1 , 6 

554  WRITE!  1 )  (  P(  HJ .  KG)  .  IIJ3  1 ,  HP  AH) 

CLOSE! UNIT3 1) 

RETURN 

EIID 


SUBROUTINE  SELF!  AF , SF , CF , AG, BG, CG, FEE) 

.'SQ=  AF$AF +BF$BF+CFSCF 
35-1=  AGS  AU+  C&.:EG+CGtfCG 
ABB=  AFSAG+BFsBC+CFsCG 
AD32= ADB+ADB 

ASAS=<  AF*BC-BF*AC)  **2+(  CF‘“-EG-BF*CG)  sx2+(  AF*CG-BF«AG)  **2 
FF=0.O0 
BO  13  HK*  1 .  10 
BO  13  NICM.HK 
LA2I»NK-MK 
A2SQ-*  ASQSLA2 l::=LA2 1 

bo  is  :il=  i ,  i i-rcc 

BO  13  HL=1,U-IIL 

lb2i=nl-ijl 

IF(  LA2 1 .  [12.  0)  GO  TO  3 
1 F  (  LB2 1 . LT . 0 )  GO  TO  3 
GO  TO  13 

3  n=SOJlT(  A2Sf'.+AD32*LA21*L321+BSQ*L321*L321) 

FF=FF+1.0/R 
13  CONTINUE 

FEE=  FF"ASAS*0 . 002 

RETURN 

END 


SUBROUTINE  MATIN!  NPAII) 

C  INVERTS  MATRIX 

DUEL'S  1011  E(  120.  120)  .  3B(  120)  ,  EST(  120) 

OPEN! UNIT*2! , FILE* ’ S’ , ACCESS* * SEQIN’ , DEVICE® ’ DSK: ’) 

DO  120  J=  I ,  IIPAN 
120  REAR! 21 ) ! E! J , I ) , Is 1 , NPAN) 

CLOSE!  UN  I T*  2 1 ) 

OPEN!  UHIT*21 , FILE*  *  E‘  ,  ACCESS*  •  SEQOUT*  , DEVICE* ’ DSIO  ’ ) 

DO  130  J*  1 .  IIPAN 
DO  11  MM* 1. NPAN 
EST< MM) *0.00 

ii  db( ran *0.00 

BB< J)= 1.0 

EST( J)* 1.0/E( J, J) 

DO  17  NIT*1.6 
DO  17  X* 1 , NPAN 
B=QB(IO 

no  15  1*1,  HP  All 

13  IFCI.IIX.X)  B=  E-E!  IC,  I )  sESTC  I ) 

EST(I0*B/'EiK.:0 
17  CONTINUE 

TRITE!  2 1 )  (  E3T!  IO  , IC*  1 .  NPAN) 

ICO  CONTINUE 

CLOSE! UN1T=21) 

TV?E  3 

3  FORMAT!’  ILYTINV  DONE’) 

RETURN 

END 

SUBROUTINE  C-Ef  IE .  YF .  ZF ,  J ,  VI ,  V2.  VO ,  VI R,  V2R.  V3R.  NBT) 

COMMON/BD^XPAN! 120) , Y?AN!  120) , ZPAN!  i 20) .AREA!  120) ,3T1  120) , 

.ACn!  120) , ACNU!  120) .AN! 120,3) , E(  120) ,?! 120.0) ,?RFS! 120) ,STOLD( 120) 
. ,PX!  120,6) 

C0NK0N/-BD2/X?T(  130)  ,  V?T!  130)  ,ZPT(  150)  .  WRF!  150)  ,  IfilFX!  150) 

.  .ICC!  150,4) 

COMMON/- ARE/RR!  SCO)  ,  XZJ!  200)  ,  YXJI200)  ,  ZYJ!  200) 

DIMENSION  X3A(2,4) ,XFA!3) ,XSAR!3,4) 

J4*J«4 
V1=0.00 
V2*  0 . 00 
V3=0 . 00 
V1R=0.0 
V2R=0.0 
VCR* 0.00 
XNJ=AN! J, 1) 

YNJ*AN! J.2) 

ZN.T=AN!  J.3) 

NS  IDE* 4 

IF!KK!  J.4)  .EQ.O)  NS  IDE*  3 
DO  20  JJ= 1, NS  IDE 
J2=  1 

IF! JJ.LT.NSIDE)  J2* JJ+  1 
J4= J4+ 1 
KF=KX<  J, JJ) 

AF=XPT! KF) 

3F=YPT( KF) 

CF=ZPT!  KF) 

R*RR!  J4) 

KG*KX! J, J2) 

ANX*  (  AF-XPT!  KG) )  /R 
ANY* ! CF-YPT!  KG) ) /R 
ANZ* ! CF-ZPT!  KG) ) /R 


TX=  XZJ(  J  >  wAHZ-YXJ(  J)  *ANY 
T/=  YXJ  <  J )  sAIIX-ZYJ  (  J ) *  ANZ 
TZ=  ZVJ (  J)*AJIY-XZJ(  J)  *AHX 
EX  1 = A* ANX+  E“  AHY+  AIIZ 
CALL  GOC  EX1 ,  R..FF,  URFC  XF)  ,  WRFC  XG)  ) 
Vl=Vl+FF*TX 
V2=V2+FF*TY 
V3=V3+FF*TZ 
XSiU  1 .  JJ)  =-.VWF(  EF) 

XSA< 2,  JJ)  =-B/UF,FC  IS') 

XSA<  3. JJ) =-C/URFC  KF) 

EX 1 n= EXl +2 . 0*ZF*ANZ 
CR--CF-ZF 

CALL  C0(  EXIF.,  R,  FH,  WRFRC  XT)  ,  I'M  KG)  ) 

VlR=Via-FE.::TX 

V2R=  V2R-FR.-TY 

V3R=V3R+FRsTZ 

IS  ARC  1 , JJ) s-A/WHFRC XF) 

X3  ARC  2 ,  J  J )  =  -  B/UXFRC  KF) 

IS  ARC  3 ,  JJ)  =CR/URFRC  KF) 

20  CONTINUE 

G= 6. 2S3 125307 

IFC  J.EQ.NBT)  GO  TO  84 

CALL  SOLID  C  ISA, G, NS  IDE) 

AGG=  A:"XNJ+B~  YNJ+C'AZNJ 
G= - S I GH  C  G , AGG ) 

84  CONTIIIUE 

CALL  SOLID  C  XSAR,  GR,  IIS  IDS) 

AGGR=  A*  XII J  D«  YII J -CR:“ZN  J 
GRsSIGII(  GR,  AGGR) 

03  CONTimJE 

7371  FORMAT!  *  G,GR» ’ , 2F15.3) 

vi=vi+xrrj*G 

V2= V2+YNJSG 
V3*V3+ZHJ=SG 

vir=vir+xnj»gr 

V2R=  V2R+ YN J”GR 
V3R=V3R-ZNJ:.“GR 

3390  FORNATC ’  VI , V2, V3= ’ ,CF13 .3) 

5591  FORMAT!  ’  V1R,V2R,V3R= ’ .3F13.5) 

RETURN 

END 
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subroutine  soliocxpn.g, "site) 
nirEITSlOH  CS<  4) , SIM  4)  ,Z(4)  ,:<PN(3.4) 

G=-6.283U)300U 

.\cni2=:a*«(  i.  i ) *mpiu  i,2)+r3’."i(2,  t)«:a>N(2.2)+xPN'.3,  n*ja>n(3,2) 

.'.c:ii3=xpr/<  i,  dsxpim  !,0)+:m(2,  i>*RPN<2.3)+xpmo,  n*xpn<3,3) 

ACR23® XPIU  1 ,2)*XPN<  1 , 3)  +'{PN(  2 , 2) *XP!« 2 , 3)  +:a5IU  3, 2)  *E?N( 3, 3) 

:F(T1S!DE.E2.4)  CO  TO  40 

C=-3.  141592659 

C3(  1>=ACR23-.ACR13:SACI112 

CS(  2) = ACR13-ACR12*ACR23 

C.S(  3)  s  ACR12-ACI123*ACR13 

sm  i)  s  jjphc  i.i)  *<  rpik  2 , 2)  *xptm  3 . 3)  -xpim  3,2)  *xpn(  2 , 3> )  + 

+  >a>ri(2,  i>*(::pn<3.2)*:iph(  i,o>-xpr«  i,2>*xpn(3,a>)+ 

+  JO’mc,  l)#CiPm  1,2)»:£PN<2.3)-XPN<2.2)*XP?U  1,3)) 

SN<2)=SN(  1) 
sn<3)=srH  i) 
sri(4)=o. 

CO  TO  30 

40  ACR14=I£PN(  1.  l)-:CTR(  1,4)+:£PN<2,  1)*XPN<  2,4)+:?N(  a,  !)  .R?N( 3,4) 

acr24=  :c?ru  i .  2> sypm  i , 4) -v:<p:i< 2 , 2) *xpri<  2, 4)  -kgtk  c , 2)  caw  3, 4> 

\cn34=:3»ti<  < ,  3)  *:i?W<  1 , 4)  +;i?:U2,  2, 4)  v/JIK  0,3)  *RPN(0,4) 

CS( 1)=ACR24-ACR14*ACIU2 

CS<  2) =  ACRl 3- ACR2C*ACR 1 2 

03 ( 3)  *ACR2”-ACR34*ACR23 

CS(  4)  =  ACRl0-ACn04  ••ACR14 

024l  =  ;CPH(2,2);Cl?II(3,4)-:a>TT(3,2)*:a,N(2,4) 

2242= EPNI  3 , 2)  sGEW  1 , 4)  -IIP.'U  1 , 2)  sEPNC  0.4) 

3243=  KHM  1.2)«;£PII(2.4)-:CPrf(2,2)«.'iPN(  1.4) 

3  io  i  *:&m  2.  i )  *;ipim  3 ,  oj  -xptu  3 . 1 )  t.-.;cPN(  2 .  o) 

Qi32=:a,ru3,  i>*xpn<  i.o-xpik  1,  i>*xpn<3,3) 

3l33=RPlH  1.  1>*RPN(2,3)-XPN<2.  lUGOW  1.3) 

S'lt  1 )  -  HPW(  1,  i>:.:524J+;*?R(2,  1)  *B242+XPII(  0.  1)«D243 

sru 2)  =-< rpiu  i  ,2)*aion-'a>n> 2,2) ::8ioa+::?iU3,2> *0133) 

STK  3)  =-(  XP!K  i ,  0)  *024 l+IlPIK  2 , 3)  ::0242+;P:K  3.0)  *3243) 

£ ;i •  4 )  -  iipru  i ,  4)  *3 10 1 +:<?m  2.4)  *b  1 32+xpn ( o ,  4)  *3 1 03 

50  CONTINUE 

D  TYPE  CC44,Sm  1)  ,C3(  1) 

D  TYPE  C244 , SN< 2) , C3<  2) 

D  TYPE  C244,SmO)  ,C3(0) 

D  TYPE  C344 . SIK  4) , C3( 4) 

G844  F0rJ!A7(  *  SR,  CS*  ’  ,  2?  13 . 8) 

scti=sn<  i )  +sn<  2)  +sn<  3)  +sru  4) 

IF(AB3(SUH) -GT.0.01)  CO  TO  25 
IF< ABSICSl 1) ) ,GT.ASS<SN< 1) > )  GO  TO  23 
IF< ABS< CS( 2) t ,GT. A3S( SN( 2) ) )  CO  TO  25 
IF(AB3(CS(G)) .GT.ABS(SH(3>>>  GO  TO  25 
1090  G=SUTI*,25 

IF(  "IS IDE.  EG.  3)  G=SUM*.  166066666667 
RETURN 

25  ST=SN(  RS IDE) 

BO  39  1=1, FS IDE 

IF( ( ABS(C3( I) ) .LT.9E-G) . AND.  (  A3S( SHI  15 ) ,LT. .9E-03) ) 

+  GO  TO  1090 

:F(5T*SR< I) .LT.O. )  GOTO  1090 
ST=Sri(  I) 

C2=CS< I )  ✓SCEYn  SIH  l)**2+CS<  I>*»2) 

G=G+AC0S(C2> 
oo  continue 

NETUTuN 

EUD 

SUBROUTINE  PREP(J) 


CCIINON/OD/NPAIM  120) ,Y?AN(  120) .ZPANC 120) , AR2A(  .20) . 3T(  120) , 

.  Ac:r<  *29)  ,  ACNV(  ICO)  .AIM  120, 3)  ,E(  120)  .Pi  12a,3)  ,?RFS<  120)  ,STOLD<  120) 
. ,?R( 120,6) 


t 


C0MK0TI/BD2/IiPT(  ISO)  ,Y?T<  130)  ,ZPT(  150)  ,  WRF(  130)  ,«  ISO) 
. ,KK(  150,4) 

corsiorr/ARE/Riu 5co) ,  xzji 200 ) ,  yxjc 200 ) , zyj<  200) 

ZYT=0. 0 
YXT®0.00 
XZT=0 . 00 
J4=  J:.“4 
JT=4 

I F C iaCf  J,4)  .EQ.O)  JT=3 
DO  20  JJ=  1 , JT 
J4=J4+1 
J2=  1 

IF(JJ.LT.JT)  J2= JJ+ 1 
:cr=EC(  J,JJ) 

~c®  kk<  j  ,  J2) 

AG»XPT(KG> 

CG=YPT(XG) 

CG=ZPT( KG) 

A?=XPT<KF) 

3F=  YPT<  ICF) 

CF-ZPT(KF) 

R® SQRT(  ( AF-AG)  **2+(  BF-3C)  *«2+(  CF-CG)  :S*2) 

irr=  af-xpaih  J) 

YT® BF-YPAIK  J) 

ZT®CF-ZPATK  J) 

AITX®  (  AF-AG)  /R 
AIJY®  (  CF-BG)  /R 
ARZ®! CF-CG) /R 
DOT®  AKX#:rr+AIlY*YT+ANZ*ZT 
jCT*:rr-DOT*AHX 
YT®  YT-DQT»AI»Y 
ZT* ZT- DOT- ARZ 
ZYT=  ZYT+ZT*  ARY-  AHZwYT 
YaT=  Yirr+  yds  anx-  ah  y«  xt 
JJZT+XTttARZ-AItXsZT 
Fall  J4)=R 

20  conti hue 

:czj(  J)=SlGi:(  AIK  J.2)  .  XZT) 

YJCJ<  J)  =  S  IGfK  AIK  J ,  3)  ,  YXT) 

ZYJ(  J)  =SIC1H  (AIK  J ,  1)  ,  ZYT) 

RETURN 

END 


